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ABSTRACT 


We  excunlne  a  method  for  solving  Liouvllle's  equation, 
consisting  of  successive  application  of  short  time  propagators 
which  are  evaluated  by  using  fast  Fourier  transforms.  The  method 
is  examined  numerically  by  computing  electronic  absorption 
spectra.  The  procedure  is  very  efficient  when  applied  to  the 
study  of  short  time  dynamics  of  systems  whose  quantum  degrees  of 
freedom  are  spatially  localized. 


INTRODUCTION 


I . 

We  propose  In  this  article  a  method  for  solving  the  equation 

ih  If  =  ±  {H^p  -  pH^}  ,  (1) 

which  describes  a  variety  of  physical  processes.  If  the  positive 
sign  is  taken  and  H  s  is  the  Hamiltonian,  Equation  (1)  is 

the  quantTira  Liouvllle  equation  (QLE)  which  gives  the  time 
evolution  of  the  non-equilibrium  density  matrix  p(t).  Taking  the 
negative  sign  and  leads  to  Heisenberg  equation  for  the 

temporal  evolution  of  an  operator  p(t).  The  form  with  different 
Hamiltonians  H^  ^  H  is  useful  for  the  computation  of  the 
electronic  absorption^  spectrum 

Z(w)  a  {l/-iT)Re  X  dt  e^“^e"‘’'^Trp(t)  .  (2) 

o 

Here 


p(t) 

-  exp{-iH^t/h)  p(o)exp{lH^t/h} 

(3) 

and 

P(0) 

with 

=  M(r)  exp{-3H^)/Z^ 

(4) 

2,  *  Tr{exp(-^H. )  ) 


3 


The  two  Hamiltonians 

H,  -  -(h^/2m)7  +V, (r)+E„  and  H,  =  -(h^/2m)vJ+V„(?) 

I  r  X  o  u  r  u 

(5) 

represent .the  nuclear  motion  in  the  "upper"  and  "lower"  electronic 

states  and  \i(r)  is  the  transition  dipole  moment  between  these 

states.  V^(R)  and  V^(R)  are  the  potential  energy  surfaces  of  the 

two  states  and  is  the  energy  gap  between  them.  The  Equation 

2-4 

(2)  is  a  rearrangement  of  the  usual  expression  that  gives  the 
absorption  spectrum  in  terms  of  the  Fourier  transform  of  the 
dipole-dipole  correlation  function.  Since  p(t)  defined  by  Eq.  (2) 
satisfies  equation  (1)  (with  the  plus  sign)  we  can  calculate 
absorption  spectra  by  solving  Eq.(l)  with  the  Initial  condition 
given  by  Eq  ( 4) . 

Thus  it  should  be  clear  that  the  development  of  an  efficient 
method  for  solving  Eq.  (1)  is  of  importance  for  non-equilibrium 
statistical  mechanics  (transport  theory  and  thermal  rates)  and  the 
spectroscopy  of  thermalized  molecules  (i.e.  prepared  in  an  oven  or 
imbedded  in  a  condensed  medium  held  at  constant  temperature) . 

We  propose  and  test  here  a  method  which  applies  to  the 
present  problem  ideas  developed  by  Fleck,  Morris  and  Feit  (FMF)® 
for  solving  Maxwell's  equations. 


II. 


THE  METHOD  OF  SOLUTION 


In  coordinate  representation  equation  (1)  can  be  written  as 

Ih  3p(r,r' ;t)/3t  »  { - (h^/2m) ( 3 ^/3r^ )  +  ( h^/2m) { 3 3r ' ^ )  + 

V^(r)  -  V^(r' ) }  p(r.r' :t)  (6) 

For  a  very  short  time  t  we  can  write 

p{r,r';t+T)  =  r .  x )  W*  (  r '  ;t  )p  (r ,  r  •  ;  t )  (7) 

where 

W^(r,T)  =  exp{(lTh^/4m)Vp}  exp{-lxV^^( r) /h}exp{  { lTh^/4m)7^} 

(8) 

5  * 

Is  the  "spilt  propagator"  used  by  FMF.  W^{r';T)  Is  obtained  from 

Eq.(8)  by  replacing  V,  (r)  with  V,(r),  with  V^,  and  1  with  -1. 

u  IT  T 

For  n  small  time  steps  we  have 

p(r,r';hT)  =  [W^( r , x ) ( r ' ,t ) ]”  p ( r , r ' ; 0 )  = 

(G(r,r' ;x)*[U^(r;x)U^(r;x)*]""^G(r,r’ :x)p(r,r':0)  (9) 

where 

G(r,r';x)  >  exp{ (lhx/4m) (7^-7^, ) >  (10) 

r  r 

and 

U^^(r,x)  =  exp{  (lhx/2m)7p)exp{-lxV^(r)/h}  .  (11) 

Uj^(r':x)*  Is  obtained  from  Eq.  (11)  by  replacing  V^(r)  with 

with  7  ,  and  1  with  -1. 

1  r  r 


We  compute  p(r,r';nT)  by  using  Eq.  (9).  The  more  difficult 

part  is  the  calculation  of  the  effect  of  the  operators  containing 

2  2 

exponentials  of  the  form  exp((iTh  /2m)7^}.  To  Illustrate  the 
procedure  we  show  in  detail  the  calculation  of 


X(r)  a  U^(r.T)f (r.r* )  . 


where  f(r,r')  is  a  known  function  and  U^(r,T)  is  given  by  Eq, 


(11)  . 


Using  Eq.  (11)  and  simple  manipulations  based  on  the 


representation  theory  we  can  write 


X(r)  =  /dk2dr^<r2|k2>  exp{ -i ( Rk2T/2m) }<k2 | r j>exp(-iTV{ r ^ ) /h) 


f(rj,r' ) 


Here  we  have  used  the  momentum  representation  in  which  the 
Laplacian  is  diagonal.  The  price  paid  for  this  easy 
diagonalization  is  that  we  must  perform  the  two  integrals  over  dr^ 
and  dk2 .  Fortunately,  since  <k2|rj>  and  <rlk2>  are  plane  waves, 
the  integrals  are  a  Fourier  transform  and  an  inverse  Fourier 
transform,  which  can  both  be  efficiently  performed  by  using  a  fast 
Fourier  transform  (FFT)  routine.® 


This  method  can  be  used  for  computing  the  effect  of  the 
operators  G  and  U.  Repeated  application  for  n  small  time  steps 
as  indicated  in  Eq.  (9),  gives  the  evolution  of  p(r,r':t)  from 


p(r,r  '  ;0)  to  p(r,r' ;nT) . 


The  number  of  operations  and  the  magnitude  of  the  errors  can 

be  analyzed  to  some  extent,  to  establish  the  conditions  under 

which  this  method  is  most  efficient.  The  error  made  by  using  the 

3 

"split  propagator"  for  a  short  time  t,  is  of  order  t  ;  if  we  were 

to  use  p(r,r';nT)  »  (U^( r , t )U^( r' ; t) • ) "p ( r , r ' ; 0) ,  like  in  the  path 
7  2 

integral  theory  ,  the  error  would  be  of  order  t  .  Furthermore,  if 

the  dependence  of  p(r,r':t)  on  r  and  r*  can  be  adequately 

described  on  a  N  point  grid,  then  each  FFT  requires,  roughly, 

Nln^N  operations  per  coordinate;  the  equations  discussed  here  have 

two  coordinates  per  degree  of  freedom.  The  number  of  grid  points 

M  is  L/1  where  L  is  the  length  over  which  p(r,r’;t)  is  localized, 

and  1  is  the  shortest  length  scale  over  which  the  p(r,r';t) 

changes.  Thus  for  each  time  step  the  number  of  operations 

2 

required  to  perform  the  necessary  FFTs  is  2(Nln2N)  The 

calculation  of  U  also  requires  N  multiplications  with 

exp{-lV^( r ) T/h} ,  one  for  each  grid  point.  The  exponentials  exp{- 

iTV(r)/h}  are  calculated  once,  at  the  beginning.  Since  two  such 

multiplications  are  required  per  time  step,  the  total  number  of 

such  operations  is  2N.  The  same  number  of  multiplications  are 

2  2 

required  for  the  terms  exp{-iTh  k^/2m) ,  where  k^  are  the  grid 
points  in  the  momentum  space;  these  exponentials  are  also  computed 
once,  at  the  beginning,  since  the  grid  points  are  not  changed 
during  the  calculation. 

The  total  number  of  operations  per  time  step  is  thus 


1  H.I  U^M  i. _l  l_.  I IPII  f> 


7 


2 

2(Nln2N)  +  4N.  The  number  of  time  steps  is  determined  by  the 
physical  characteristics  of  the  problem  at  hand.  In  the  case  of 
electronic  absorption  spectrum,  which  is  used  for  Illustration  in 
this  article,  the  time  step  must  be  smaller  than  2it/6,  where  6  is 
the  width  of  the  Franck-Condon  envelope,  and  the  length  of  time 
that  p{t)  must  be  propagated  is  of  order  y~^,  where  y  is  the  width 

O 

of  the  narrowest  line  that  we  want  to  resolve. 


III.  The  Calculation  of  Electronic  Absorption  Spectra 

To  test  the  method  described  here  we  have  applied  it  to  the 
problem  of  electronic  absorption  by  a  model  diatomic  molecule 
whose  ground  and  excited  state  potential  energies  are  harmonic 
oscillators  of  equal  frequency  =  0.01  a.u.  =  2.72  meV  =  2176 
cm  whose  equilibrium  positions  are  displaced  by  a  length  6x  = 
0.33  a.u.  ^  0.16  A,  and  whose  reduced  mass  is  that  of  a  proton. 

9 

We  chose  this  problem  because  the  exact  result  is  known. 
Since  the  spectrum  and  the  Hamiltonians  have  all  the  features 
expected  in  a  realistic  problem,  the  model  provides  an  adequate 
test  of  the  method  of  computation. 

The  calculation  starts  with  the  assumption  that  the 
transition  dipole  ^(r)  is  independent  of  r.  There  is  no 
additional  difficulty  in  using  a  function  for  ^(^);  we  used  a 
constant  since  exact  results  are  available  for  that  case.  The 
equilibrium  density  matrix.  exp{-3Hj^)/Zj^  was  calculated  in 
coordinate  representation  by  using  a  FPT  method  proposed  by 
Hellsing,  Nitzan  and  Metlu.^^  This  provides  the  exact  initial 
value  p{r,r';t=0),  »  <r  |  exp( -pH^^)  |  r  '  >/Z  ^  .  This  is  propagated  by 
using  (9)  and  the  method  described  above.  The  spectrum  Z(u)  is 
obtained  from  Eq.  (2).  Since  the  band  edge  is  determined  by  the 
energy  gap  E^  between  the  states  (i.e.  the  electronic  excitation 
energy)  we  can  shift  the  spectrum  arbitarily  on  the  frequency 
scale  without  causing  confusion. 

The  length  of  time  T  for  which  we  need  to  propagate 
p(r,r’;t)  is  set  by  the  magnitude  of  r  appearing  in  Eq.  (2).  Her 


y  is  the  rate  of  the  slowest  process  causing  the  decay  of  the 
upper  state  amplitude,  and  it  could  be  the  natural  line  width, 
the  rate  of  a  radiationless  transition,  predissociation,  or  photo¬ 
dissociation  in  an  additional  degree  of  freedom  whose  dynamics  is 

4  12 

not  explicitly  Included;  '  one  can  also  think  of  exp(-7't)  as  a 
"filter  function"®  used  to  avoid  the  numerical  difficulties 
associated  with  computing  6  functions  numerically.  In  the  latter 
case  y  must  be  smaller  than  the  transition  frequencies  which  one 
intends  to  resolve  when  the  spectrum  is  computed.  We  use 

throughout  the  paper  yT  =  5  with  y  =  0.028  a.u.  (y/Wq  =  0.28). 

The  dependence  of  the  spectrtim  on  the  number  of  spatial  grid 

points  is  illustrated  in  Fig.  1.  The  result  obtained  with  32  grid 

points  is  exact,  that  for  16  points  is  reasonably  accurate,  while 
the  use  of  eight  points  leads  to  serious  errors. 

The  sensitivity  to  the  number  of  spatial  grid  points  appears 
because  the  height  of  each  peak  and  the  total  width  of  the 
spectrum  is  determined  by  the  magnitude  of  the  Frank-Condon 
factors,  which  are  wave  function  overlaps,  A  poor  spatial  grid 
will  result  in  erroneous  high  energy  Pranck-Condon  factors,  since 
the  highly  excited  wave  functions  appearing  in  them  have  a  larger 
number  of  spatial  oscillations  and  these  are  not  well  described  on 
a  coarse  grid.  When  an  insufficient  number  of  grid  points  is  used 
the  low  frequency  part  of  the  spectrum  is  usually  not  as  bad  as 
the  high  frequency  one. 

In  Figure  2  we  show  the  dependence  of  the  spectrum  as  a 
function  of  the  number  of  time  steps.  The  results  obtained  by 


using  1000  time  steps  are  exact,  those  using  25  steps  are  adequate 
and  the  error  made  by  this  dramatic  lowering  of  the  number  of 
steps  occurs  mainly  at  the  low  intensity  side  of  the  spectrum, 
where  the  experiments  are  also  likely  to  be  less  accurate. 

In  Figure  3  we  show  the  temperature  dependence  of  the 
spectra.  The  Increase  in  temperature  causes  the  growth  of  the  hot 
bands  and  a  corresponding  decrease  in  the  peaks  which  are  present 
at  zero  temperature. 

IV  DISCUSSION 

Both  actual  calculations  and  estimates  of  the  number  of 
operations  indicate  that  the  present  method  of  solving  the 
Liouvllle  equation  is  very  efficienct  if:  (a)  the  density  matrix, 
in  coordinate  representation,  is  spatially  localized;  (b)  the 
density  matrix  changes  with  position  over  a  sizable  space  scale; 

(c)  the  duration  of  the  dynamic  process  of  interest  is  short  (in 
spectroscopy  this  means  that  the  Pranck-Condon  envelope  is  broad) ; 

(d)  the  required  time  step  is  small  (in  spectroscopy  this  means 
that  the  ratio  between  the  width  of  the  narrowest  line  and  the 
width  of  the  envelope  is  large);  (e)  the  problem  has  few  quantum 
degrees  of  freedom. 

Even  though  a  theory  based  on  Liouvi lie's  equation  is  more 
general  than  one  provided  by  the  Schrodlnger  equation,  such  a 
formulation  is  not  always  required,  nor  is  it  generally  advisable. 
In  cases  when  k^T  is  comparble  to  the  excitation  energy  of  the 
quantum  system,  it  is  computationally  more  efficient  to  propagate 
the  few  states  which  are  thermally  populated  and  average  the 


results  with  a  quantum  Boltzman  factor.  This  happens  because  the 
dimensionality  of  the  Liouville  equation  is  twice  that  of 
Schrodinger  equation. 

When  we  deal  however  with  molecules  imbedded  in  a  condensed 
medium,  whose  molecules  are  held  at  constant  temperature  and  are 
propagated  by  classical  mechanics,  the  situation  is  not  so  clear 
cut . 

Under  these  circumstances  the  initial  conditions  for  the 
quantum  degrees  of  freedom  are  determined  in  two  steps:  First  we 
use  the  Monte  Carlo  method  to  establish  the  position  of  the  atoms 
of  the  medium  and  to  compute  the  potential  energy  of  molecule  in 
the  presence  of  the  medium;  then  we  determine  the  quantum  density 
matrix  of  the  molecule  for  these  potentials.  Under  these 
conditions  the  molecular  wave  functions  are  much  more  difficult  to 
determine  than  the  equilibrium  density  matrix.  Therefore,  since 
the  initial  state  is  defined  by  the  density  matrix  the  subsequent 
time  evolution  (caused  by  the  electromagnetic  field  and  the 
classical  motion  of  the  molecules  of  the  medi\im)  and  the  spectrum 
corresponding  to  it,  must  be  calculated  by  solving  the  Liouville 
equation.  The  same  is  true  for  rate  or  transport  processes. 

The  reason  why  the  present  method  might  play  an  important 
role  in  future  developments  is  that  it  is  so  efficient  that  such 
calculations  -  which  require  the  calculation  of  the  quantum 
^volution  of  the  molecule's  density  matrix  for  each  Monte  Carlo 
configuration  of  the  molecules  in  the  medium  -  are  possible  if  the 
problem  has  one  or  perhaps  two  degrees  of  freedom  which  are 


spatially  localized,  and  the  process  is  completed  In  a  reasonably 
short  tine. 
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FIGURE  CAPTIONS 


Fig.  1. 


Fig.  2. 


The  dependence  of  the  computed  spectrum  on  the  number 

of  points  N  in  the  spatial  grid.  Full  line  N>8 ; 

dotted  lines  N«16;  dashed  line  N«32.  We  have  used 

1000  time  steps,  *  12.5,  «  »  272  meV  *  0.01 

Op  o 

a.u. ,  y/w^  =  0.28. 
o 

Spectra  computed  with  various  time  steps:  full  line, 
1000  steps;  dashed  line,  100  steps;  and  dotted  line, 
25  steps.  The  spatial  grid  has  64  points.  The 
parameters  are  those  used  in  Figure  1. 


Fig.  3. 


Spectra  at  various  temperatures.  The  parameters  are 
those  used  in  Fig.  1. 
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